Spacetime thermodynamics and subsystem observables in a kinetically constrained 

model of glassy materials 
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In a recent article [M. MeroUe et al., Proc. Natl. Acad. Sci. USA 102, 10837 (2005)], it was argued 
that dynamic heterogeneity in d- dimensional glass formers is a manifestation of an order-disorder 
phenomenon in the d+1 dimensions of space-time. By considering a dynamical analogue of the free 
energy, evidence was found for phase coexistence between active and inactive regions of space-time, 
and it was suggested that this phenomenon underlies the glass transition. Here we develop these 
ideas further by investigating in detail the one-dimensional Fredrickson-Andersen (FA) model in 
which the active and inactive phases originate in the reducibility of the dynamics. We illustrate the 
phase coexistence by considering the distributions of mesoscopic space-time observables. We show 
how the analogy with phase coexistence can be strengthened by breaking microscopic reversibility 
in the FA model, leading to a non-equilibrium theory in the directed percolation universality class. 
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I. INTRODUCTION 

The hallmarks of the glass transition j| are a sudden 
increase of relaxation times, and dynamic heterogeneity 
. A recent paper Q described the idea that these phe- 
nomena are manifestations of phase coexistence in space- 
time. The two competing phases are an active state, 
containing many relaxation events, and an inactive one, 
where these events are scarce. The purpose of the cur- 
rent paper is to develop these ideas by building on the 
results of Ref . [3j] . 

In the active phase of the dynamics, the existence of 
the inactive phase can be inferred by measuring the dis- 
tribution of any observable that measures the quantity of 
relaxation events within a finite space-time window [3|. 
Because of the coexisting inactive phase, these distribu- 
tions display non-Gaussian tails. This behaviour is anal- 
ogous to the situation in standard phase transitions: the 
distribution of cavity sizes in a fluid near liquid vapour 
coexistence has non-Gaussian tails [4], as does the distri- 
bution of box magnetization for the Ising model in the 
vicinity of its phase transition [5| . 

The concept of space-time thermodynamics that we 
use Q is illustrated in Fig. [1] It employs a simple pic- 
ture of glassy dynamics at low temperatures, in which 
diffusing excitations coalesce and branch Q. This dy- 
namics is reducible [7J]. That is, there are two steady 
states, an active one with a finite density of excitations, 
and an inactive one with strictly zero density. Working 
in the active state, the dominant fluctuations on large 
length scales are space-time regions in which there are 
no excitations, as illustrated in Fig. [TJ We refer to these 
regions as "bubbles" of the inactive phase. 

We characterise the bubbles by the probability density 
function -Pbubbic(^ ''') for their spatial and temporal ex- 
tents (denoted / and r, respectively). To the extent that 
these rare fluctuations are dynamical analogues of those 
near phase coexistence, and for bubbles which are large 
compared with the bulk correlation length of the active 
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FIG. 1: Illustration of a trajectory in a facilitated model with 
space-time "bubbles" of the inactive state. The boxes illus- 
trate finite observation space-time windows: the top one cor- 
responds to a typical region; the bottom one is a rare collec- 
tive fluctuation of size much larger than those typical of the 
active state. On the right are trajectories from the one-spin 
facilitated one-dimensional Fredrickson-Andersen model [1G|. 
at r = 1 for observation windows of L = 180 and tobs = 320 
(smaller observation windows of L = 60 are also outlined). 



state, we can write: 

^bubbio(^,T) ex exp {- [aj + a2T + ^iIt + A{1 - vrf] } , 

(1) 
where cti and CT2 are, in effect, surface tensions in the spa- 
tial and temporal directions respectively, /z is the differ- 
ence between the bulk free energy densities of active and 
inactive phases, and A controls the aspect ratio of the 
bubbles, in conjunction with the velocity parameter v. 
This distribution of ^ can be investigated by consider- 
ing the probability distribution functions for observables 
that are averaged over finite regions of space-time, such 
as the boxes in Fig. [T] [J] . Large rare bubbles dominate 




FIG. 2: Sketch illustrating the effect of the choice of initial 
condition in a model of diffusing excitations that branch and 
coalesce. An initial state with no excitations (left) persists 
throughout the observation time. All other initial conditions 
result in the system exploring the active steady state (right). 



the tails of these distributions. 

The terminology "free energy" , as used in the previous 
paragraph, refers to space-time; the relevant ensemble 
is a the space of histories (or trajectories) of the system, 
which is observed for a specified period of time, iobs- This 
concept of a free energy is central to Ref. [3|, and it is 
found also in the thermodynamic formulation of dynam- 
ics due to Ruelle p, !3|. In particular, for a d-dimensional 
system, one considers statistical mechanics in the corre- 
sponding [d + l)-dimensional space. The extra dimen- 
sion is that of time, so the resulting space is necessarily 
anisotropic. 

Thus, to appreciate the analogy with phase coexistence 
for d = 1, consider an anisotropic two-dimensional Ising 
model below its critical temperature, with a small field 
stabilising the spin-up phase. The spins are predomi- 
nately up, but there will be a distribution of spin-down 
domains with various sizes. In the case where the inter- 
actions are different in the x and y directions, the distri- 
bution of sizes for domains of the spin-down phase is of 
the form 

-Pdown(a;, y) oc exp {- [a^x + cTyy + ^i'xy + K'{x - v'yf] } 

(2) 
where x and y are the extents of the bubbles in the x 
and y directions; a^ and <jy are surface tensions; /z' is 
the free energy difference between up and down phases 
(proportional to the small field); and A' and v' set the 
shape of the bubbles. 

As phase coexistence is approached in the Ising sys- 
tem, the free energy difference /i' approaches zero. In 
the dynamical system, the corresponding limit is /i — *■ 0, 
in which case active and inactive states coexist. A cen- 
tral result of this article is that /i = for systems in 
which the dynamics is reducible into active and inactive 
partitions. Within the thermodynamic formalism of Ru- 
elle the free energy difference /i is the difference in the 
so-called 'topological pressure' [1, [0] between the active 
and inactive phases of the dynamics. 

Since the phases are at coexistence in these models, a 
boundary condition suffices to make the system choose 
one phase. For example, consider a system with periodic 
spatial boundary conditions, as sketched in Fig. [51 If 
we choose our initial condition to be completely inactive. 



then this state persists forever; all other initial conditions 
lead to the active state. 

In this article, we use the the Fredrickson-Andersen 
(FA) facilitated spin model [10] as an illustrative system. 
In the following section, we probe the distribution of ([1]) 
by analysing the distributions of space-time observables, 
such as the total amount of dynamical activity inside a 
given spatial region, integrated over a finite observation 
time. The forms of these distributions also imply that 
/i = 0, in accordance with our analytic arguments. 

We also compare the FA model with two other sys- 
tems. In Ref. 11], a model with appearing and annihi- 
lating excitations (the so-called AA model) was shown 
to have the same two point correlations as the FA model 
(up to a multiplicative factor). However, we show that 
this model has mesoscopic fluctuations that are very dif- 
ferent from those of the FA model. Conversely, breaking 
microscopic reversibility in the FA model moves moves 
the system into the directed percolation (DP) universal- 
ity class [13 • In that case, the the two point functions 
and critical scaling change qualitatively. However, the 
active and inactive states still coexist, so the fluctuations 
are still well-described by the phase coexistence analogy. 

The paper is organized as follows. In section |TT] we de- 
fine the model, discuss its trajectories, and introduce the 
observables of interest. We also show that reducibility 
of the dynamics is essential to see the phase coexistence 
effect. In section iHll we formalise our discussion of space- 
time thermodynamics by considering the distribution of 
the dynamical action that is analogous to the thermo- 
dynamic free energy. We consider the effect of breaking 
microscopic reversibility in section [TVl 



II. MODELS, TRAJECTORIES AND 
OBSERVABLES 



Kinetically constrained models (see Ref. [7| for a 
comprehensive review) are defined in such a way that 
their non-trivial behaviour is purely dynamical in ori- 
gin. Thus, they are the natural framework for inves- 
tigating dynamical heterogeneity and its consequences 

[aEtilliiEi,^- 

A very simple kinetically constrained model is the 
single-spin facilitated FA model [3, [l3| in dimension 
d = 1. It is defined for a chain of N binary variables 
n.i G {0, 1}, with trivial Hamiltonian H — J^^Ui. The 
model evolves under the following Monte-Carlo (MC) dy- 
namics. At each MC iteration a randomly chosen site i 
changes state according to: 






rii = 1 probability /^e ^ ^ 
rii — Q probability /i. 



where we have set J — 1 and (3 = T ^. The non-trivial 
part of the dynamics is due to the facilitation function 



fi = n-i+i + "i-i - ?^^+l?^^-l• 



(3) 



That is, a spin flip on site i can take place only if at 
least one if its nearest neighbours is in the excited state. 
Since fi does not depend on Ui then the time evolution 
of the model obeys detailed balance with respect to H at 
temperature f3~^. The mean density of up spins in the 
active state is 



€={1 + 6") 



(4) 



The unit of time is an MC sweep {N attempted spin 
flips). 

In the following, we consider an FA model in which 
N is to be taken to infinity in the thermodynamic limit. 
We will take a subsystem of this model, containing L + 2 
spins, {no, . . . ,ni_|_i}. The length L is to remain finite: 
it is a mesoscopic quantity. A trajectory for the subsys- 
tem specifies the state of the {L + 2) spins at iViobs MC 
steps. The trajectory is defined within a space-time ob- 
servation box of size Ltohs- We consider observables such 
as the density of excitations within the box, for a given 
trajectory: 



TOtraj = (LNtohs) ^ ^ 



= 1 T = l 



(5) 



where riir is the state of the ith spin in the rth state of 
the trajectory. Note that the boundaries uq and ri^+i do 
not appear in the sum. 

The quantity TTitraj is analogous to the net magneti- 
zation of a (d + l)-dimensional Ising model, and thus 
we often refer to it as the "magnetization". However, 
its physical meaning is that of mobility or activity, i.e., 
the net number of cells in the observed space-time vol- 
ume that exhibit an appreciable degree of molecular mo- 
tion. This meaning is that of Garrahan and Chandler's 
coarse grained view of structural glasses [g, |Tj|. In 
an atomistic treatment, a corresponding quantity can 
be expressed in terms of the observable Fj{k,dt,t) ~ 
cos{k-[rj{t+St) — fj{t)]} where 'rj(i) is the position of the 
jth particle at time i, and fc is a microscopic wave vector 
(on the order of the inverse particle size). If Fj{k,St,t) 
is small, then particle j experienced a relaxation event 
between times t and t + St. For a subsystem of parti- 
cles, 5, we define the density of relaxation events, or net 
mobility, in space and time to be 



Ktr^iik,6t)^ (tobsNs) 






dt[l- Fj{k,St,t)] 



(6) 

where Ns is the number of particles in the subsystem S. 
The quantity i^traj(fc, St) is parametrised by microscopic 
length and time scales (|fc|~^ and St), and is summed 
over a mesoscopic region of space-time. The parameters 
|fc|~"'^ and St coincide with coarse graining lengths and 
time, respectively. 

The fluctuations in KtmjikjSt) are related to four 
point observables |T6[; the distribution of KtTa.j{k,St) is 
also related to the distributions of correlation functions 
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FIG. 3: Distribution of trajectory magnetization P(m) at 
/3 = 1, 1/ = 60 and various observation times. We use A'^ = 180 
which is large enough so that P(m) does not depend on A'^. 
The exponential tails of P(m) all have similar gradients: the 
dotted lines are P{m) oc exp(o"m) with a = 87. 



measured in [T^, [I^ . Results for distributions of the den- 
sity of relaxation events in atomistic glass-formers show 
non-Gaussian tails similar to those that we present here 
[20|. However, the purpose of this article is to discuss 
generic properties of dynamically heterogeneous systems, 
so we restrict ourselves to very simple models. 



A. Distribution of the event densities 

We consider the distribution of the trajectory activity, 
or box magnetisation rritraj '■ 



P{m) = ^ Pua^jS{m - mtraj)- 



(7) 



traj 



The sum is over all possible trajectories of the system, 
with their associated probabilities Ptraj- In the limit of 
L — > cxD or tobs — *■ oOj the central limit theorem states 
that P{m) will be sharply peaked around m = c, with 
a variance that scales as (Ltobs)"^- The observation of 
Ref. @ was that while this is the case in the limit of 
large (Ltobs), the distribution of m remains non- trivial 
even for L^ ^ and iobs ^ t, where ^ « c~^ and r « c~^ 
are the correlation length and time associated with the 
zero temperature dynamical fixed point of this model. 

We show various P{m) in Fig. [3l The peak of P{m) 
is Gaussian, and this corresponds to trajectories like the 
one on the top-right of Fig. [1] The tails of P{rn) are 
exponential for small m. This data was obtained using a 
combination of transition path sampling (TPS) [2l| and 
umbrella sampling [22]. For details see the Appendix. 
The trajectories in the tail are dominated by rare large 
regions which have no excitations, like the the one on the 
bottom-right of Fig. [1] 

To explain the presence of the exponential tail, sup- 
pose that we have a trajectory with an initial condi- 
tion containing a large empty region of size a;, and that 



this empty region persists throughout iobs, leading to a 
box magnetisation that is small. Now consider a tra- 
jectory whose initial condition has an empty region of 
size X + (Aa:), but is otherwise identical to the original 
one. Then the probabilities of the two trajectories are 
in the ratio (1 — c)'^^, and their magnetisations differ by 
approximately c{Ax)/L. Thus we predict 



P{m) - (l"c)' 



ih I C 



exp[(mL/c)ln(l-c)]. (8) 



The significance of this result is that the right hand side 
is independent of iobs^ increasing the observation time 
changes the factor multiplying the tail in P(m\ but the 
gradient of the tail remains constant. Fig. [3] shows that 
the gradient is indeed independent of iobs- However, the 
quantitative prediction for its value is accurate only to 
within 20-30%. This inaccuracy arises because the large 
inactive region may extend beyond the edge of the box. 
Thus, adding extra down sites to the initial condition 
may decrease the magnetisation by an amount less than 
{cjL). 

In Fig. [3l the exponential tail of Pirn) depends only 
on L and not on tobsj but this is a result of the shape of 
the observation box: the behaviour of P(rn) is symmetric 
with respect to L and tobs- In Fig. [3] we show how P{m) 
changes as we change the aspect ratio of the observation 
box. There is a crossover at L ~ wiobs, where v ~ 0.06 
a,t (3 — 1. For L > wiobs the gradient of the exponen- 
tial tail is independent of tobs as described above, but for 
L < vtohs then the argument must be modified. Typical 
trajectories at small m change their form to that shown 
in Fig. 2] (right). In that case, increasing the width of the 
bubble does not change m; we must instead increase the 
time for which the inactive state persists. Repeating the 
argument leads to a situation in which the gradient of the 
exponential tail is proportional to iobs and independent 
of L, see Fig. 2] In other words, the gradient of the expo- 
nential tail is proportional to max(L,tiiobs)- The spatial 
and temporal extent of the box enter on equal footing. 
The data of Fig. H] is consistent with the crossover from 
large L to large iobs occuring when tobs is of the order of 
the first passage time across a bubble of size L. This time 
scales as (L/v) ~ {L/c^) since the spreading velocity [I4I 
in the FA model is w ^ (^/''") ^ c^- 

Consistent with the known scaling of the FA model 
in one dimension [7], we find that the temperature de- 
pendent of P{m) can be accounted for by using rescaled 
variables. The function 
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FIG. 4: We show P(m) at /3 = 1 for varying L and tobs. 
We use A'' = max(120,3I/) which is large enough that the 
results do not depend on N . (Top) Increasing observation 
time at at fixed L — 20. (Middle) Increasing box size L at 
tobs = 320. As L or tobs is increased, we move from a regime 
in which the tail gradient is independent of the increasing 
parameter to a regime in which the gradient is proportional 
to that parameter. (Bottom) We show a typical trajectory for 
large tobs and L = 20 where the observation box is outlined: 
the size of the total spatial region shown is A'^ = 60. For large 
L the trajectories are of the form shown in Fig. [T] 



P{m/c;cL,chobs) 



(9) 



depends only weakly on c, as shown in Fig. \E\ While 
subleading corrections to the scaling do appear to be sig- 
nificant, there is no qualitative change on lowering the 
temperature. This justifies our working at inverse tem- 
peratures around /? = 1, where the data is easier to obtain 
than in the more glassy regime of the FA model (/3 > 1). 
As a final observation in this subsection, we note that 
the behaviour of P{m) for small observation boxes, cL < 



1 or c'^tobs < 1, is different from that for large boxes. This 
is shown in Fig. [G] The distribution acquires a secondary 
peak at m — > 0: there is a finite probability of observing 
TO = 0, so that P(m) diverges as ttt- — > 0. We will return 
to this feature below. 



B. Comparison with a model of appearing and 
annihilating excitations (AA model) 

In the previous subsection, the non-trivial structure 
in P(m) occurs in the tails of the distribution. It is 
associated with rare collective fluctuations in the sys- 
tem. We will argue below that these rare fluctuations 
do contain significant information about the model. We 
first compare the magnetisation distributions of the FA 
model and the AA model of [11| . The AA model has the 
same two-point correlation functions as the FA model in 
equilibrium, but other properties of the two systems are 
qualitatively different. 

We define the AA model by specifying local moves for 
a chain of binary variables: 



{ni,ni+i) = (1,0) 
{ni,ni+i) = (1,1) 
{ni,ni+i) = (0,0) 



(nj,ni+i) = (0, 1), rate 76 

{nt,ni+i) == (0,0), rate 7 

(nj,ni+i) = (1,1), rate 76' 



-0' 



-213' 



where 7 = 2/(1 — e~^ ) is an arbitrary temperature 
dependent factor that rescales the time, chosen for later 
convenience. These rates again respect detailed balance 
with respect to H = ^^ n^, and the inverse temperature 
is/3'. 

Correlation functions in the steady state of the AA 
model at a given temperature are related to those in the 
FA model at a different temperature [ll| (strictly this 
holds when the facilitation function is proportional to the 
number of up neighbours, not when it is equal to zero or 
one as here, but there is no qualitative change to the 
physics). The relationship between correlators depends 
only on simple multiplicative factors, for example ^11,] 



{nitnjt')FA,i3 
4 



)fA,/3 



1 



—3 i{nttnjt')AA,i3' - (n-)AA,/3') 



(10) 



where /3 is the inverse temperature in the FA model and 
/3' the inverse temperature in the AA model. Equation 
(fTU)l holds when 



Vl + e-/3 - 1 



Vl + e-/3 + l 



(11) 



It is clear from ([T0| that the scaling behaviour and crit- 
ical properties of the FA and AA models are the same. 
However, their activity or magnetisation distributions are 
different, as shown in Fig. [T) The exponential tail is ab- 
sent in the AA case. Moreover, the persistence functions 
of the two models are different [the persistence function is 
the probability that a site does not flip at all in an inter- 
val of length t] . The fluctuations responsible for the tails 
in P{m) can be linked to the decoupling of exchange and 
persistence times in the FA model [23], a feature which 
is absent in the AA case. We argue that there are im- 
portant differences between the dynamics of the FA and 
AA models: subsystem observables like P{m) make these 
differences clear where two-point functions like (riiT^-jo) 
do not. 
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FIG. 5: Data showing (approximate) scaling of P{m) in the 
FA model at various temperatures, scaled according to ([S}. 
We plot P{m/c): the box sizes are L = (60,88,135); the 
observation times are tobs = (320,1280,5330); and we use 
A^ — 31/. These temperatures are not very small, so there are 
subleading corrections to scaling, but there is no qualitative 
change to the scaled distribution on lowering the tempera- 
ture. Further, the computational time required at /3 = 2 is 
quite significant, so we cannot rule out small systematic er- 
rors arising from non-convergence of our TPS procedure (see 
appendix) . 
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FIG. 6: Plot of P{m) at /3 = 1, L = 2, iobs = 160, showing 
secondary maximum at small m. (Inset) Enlargement of the 
secondary peak, shown in a linear scale for P(m). 



III. DYNAMICAL ACTION AND 
THERMODYNAMIC ANALOGY 



We have shown that the distribution of subsystem 
magnetisations in the FA model has an exponential tail at 
small magnetisation. We also showed that this tail comes 
from rare mesoscopic regions of space-time in which there 
are no relaxation events. In this section we study the sta- 
tistical mechanics of the space-time configurations of the 
subsystem. We draw an analogy between these space- 
time statistics and the thermodynamic statistics of a sys- 
tem near a phase coexistence boundary. 
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FIG. 7: Distribution of (reduced) box magnetization in the 
FA and AA models. The reduced variable x — {m — {m))/am 
where a^ = {{jn — (m))^) — c(l — c)L~^ is the variance of 
the instantaneous magnetisation. Parameters are L — 60, 
tobs = 320, A'' = 180; in the FA model /3 = 1; in the AA 
model /3' is given by ([TT} with /3 = 1. For the AA model, 
P(m) is close to Gaussian. The standard deviation am is not 
trivially related to the variance of the box magnetisation, so 
the fact that the Gaussian parts of the two distributions are 
very similar is a non-trivial consequence of the exact mapping 
between the two models. 



A. Thermodynamic analogy 

To make this analogy concrete, we define the probabil- 
ity of a trajectory being generated by some (Markovian) 
dynamical rules: 

^traj = W W{Sr+l\{Sr,no^r,nL + l,r}) 



T = l 



xPo(si)-Pbc("-0,t,«L + l,t) 



(12) 



where Sr = {?iir, ■ • • , nLr} represents the state of the sys- 
tem at time r; the transition probabilities are denoted by 
W{s'\{s,nQ,nL+i}); Phcino,t,nL+i,t) is the probability 
of the trajectory of the two boundary spins; and Po{si) 
is the probability of the initial condition (si). We note 
that the transition probabilities into a state St+i depend 
on the state of the boundary spins as well as the state 

Sr- 

It is convenient to define the dynamical action, 



(13) 

(14) 

which is equal to unity when b = 1. The extensivity 
properties of the action £ will be discussed below. To 
make the correspondence with a thermodynamic parti- 
tion sum, we simply define 



^traj — — In Ptiaj , 

and the dynamical partition sum, 

traj 



E' 



-/3(7co„ 



(15) 



where the sum is over configurations of a system; f/conf is 
the energy of a configuration; and /3 is the inverse tem- 
perature. The application of thermodynamic formalism 
to dynamical objects such as Zi, is found in the work of 
Ruelle [Sj] . A very recent and detailed discussion of its ap- 
plication to Markov chains is given in 9] ; the parameter 
s of that paper is 1 — 6 in our notation; the Kolmogorov- 
Sinai entropy h^g of the Ruelle formalism is the mean 
action, {£) in our notation. 

To treat phase coexistence in ordinary thermodynam- 
ics, one may decompose the partition function into its 
pure states [24 1 



Z 



Y.Z. 



(16) 



where Za is the contribution of pure state a to the parti- 
tion sum. For example, in the Ising system then the two 
pure states have positive and negative magnetisations, 
and statistical weights Z| and Z|. The difference in free 
energy density between the two states is 



m'= lim N-Hog{Z^/Z^). 



N 



(17) 



At phase coexistence, we have /i' = 0. The result is that 
subextensive contributions to log(Z|/Zj^) from bound- 
ary conditions are sufficient to select one phase over the 
other. 

By analogy, for dynamically heterogeneous systems, we 
make an decomposition of Z into its reducible partitions 
[7| . If there are two such partitions then we write 



Zb=i-Pic(l) + Pic(2) 



(18) 



where the two terms are the probabilities that an initial 
condition lies in one partition or the other. We then 
define 



M 



lim 



(iViobs)-Mog[Pic(l)/Pic(2)], (19) 



by analogy with the static free energy difference /i'. 

The probabilities Pic(l) and Pic (2) are independent of 
tobs, so PL vanishes as iobs is taken to oo, and an initial 
condition on the system is sufficient to select the active or 
inactive partition of the dynamics: see figure[2l In a sim- 
ilar way, a boundary condition is sufficient to select the 
up or down pure state in an Ising system below the crit- 
ical temperature. Recalling ([1]) and ^, we observe that 
phase coexistence and reducibility lead to distributions 
with /i' = and /i = respectively. The free energy cost 
associated with large bubbles (or down domains) comes 
from their boundaries and not from the bulk. 



B. Discussion of magnetisation distributions 

Fig. [1] illustrates how FA model trajectories with 
low densities of relaxation events typically contain a 
large compact "bubble" in which there are no relaxation 



events. In the FA model then these large bubbles are 
more common than homogeneous reductions in the den- 
sity. Conversely, in the AA model, the large bubbles are 
absent and the homogeneous regions of low density de- 
termine the tail of the magnetisation distribution (recall 
figure [7|) . 

The FA model is reducible, so the distribution of large 
bubbles will be that of (P), with fj. = [recall HH)]. 
The derivation of the magnetisation distribution from 
the bubble distribution depends on the ratio of L and 
iobs- For observation windows such that wiobs < L, tra- 
jectories with small mtraj contain large bubbles that span 
the temporal extent of the observation window, as in the 
lower right panel of figure [TJ In that case, we have 



m~c[l-(Vi)]. 



(20) 



C. Distribution of the dynamical action 

In the analogy between static and dynamic partition 
sums, the dynamical action of Eq. (|13p corresponds to 
the energy of the static system. For the FA model, the 
transition probabilities are 

1 ^ 

W{Sr+l\{Sr,nor,nL+l,T}) = Wq + Vs^+^ ,Sr J^'^Wi 

i=\ 

(23) 
where the projector Vs.s' takes the value of unity if s and 
s' differ in the state of exactly one spin, otherwise it is 
equal to zero; 

Wi = jir\{\ - 7ii,r+i)".i,T + "■j,T+i(l - nir)e~^\ (24) 
accounts for transitions into state s; and 



where I is the spatial extent of the bubble. The probabil- 
ity of the box containing a bubble of this size and shape 
is 



p(0 



drP(/,T). 



(21) 



Assuming that A is not too small, so that P{1,t) is 
sharply peaked around I = vt, then 



In P(m) 

am 



[L/c) [ai + [a^/v] -f 2{^i/v)L{\ - m/c)] , 

(22) 

for 777, < c[l — {vto\^sl L)]. Unless /i is very small, the term 
proportional to /i will dominate at large L, leading to a 
Gaussian distribution. However, in the case of small /x, 
the distribution is exponential in tti, and the gradient of 
this exponential tail is proportional to L and independent 
of tobs, as observed in Fig. [31 A similar argument holds in 
the opposite regime of L < wiobs and 777 < c[l~{L / vtohs)], 
which explains the exponential tails of Fig. [H 

Pursuing the analogy with phase coexistence, we iden- 
tify the active partition of the dynamics with a dense 
phase and the inactive partition with a sparse phase. The 
density of dynamical activity in the active phase is (777) — 
c. The typical length scale associated with bubbles of the 
sparse phase is crj~ which scales as c""'^. Thus, the typical 
bubble size and the typical spacing between excitations 
scale in the same way as the temperature is reduced. This 
scaling [ll| is consequence of detailed balance in the FA 
model, which implies that excitations are uncorrelated at 
equal times: {{riiT — c){njr — c)) = c(l — c)Sij. 

In section IIVI we will show how generalising the FA 
model to a system which does not obey detailed balance 
leads to a situation in which the typical bubble sizes are 
larger than the excitation spacing. This situation is the 
usual one in systems at phase coexistence: the FA model 
is a special case because of the extra symmetry of detailed 
balance. 



Wo = 6s^^„s^ ] 1 - ^"' II [/-(I - "-)e"^ + n^f^r] \ 

(25) 
accounts for transitions from state s. The symbol Ss^s' is 
equal to unity if and only if states s and s' are identical, 
and the function fir was defined in ([3]) : its value is unity 
if spin i is free to flip; otherwise it is zero. These operators 
enforce the constraint that only one spin flips on each 
time step. 

Taking the limit of large TV, and using (J13p . we arrive 
at the action for an allowed trajectory: 

41! = A/'upflipsln(iVe'3)+AAdownflipsln(iV) 



-N' 






f^A{l - n.rje-'' + ' 



(boundary terms). 



(26) 



where Wupflips = Z^ir "«,-r+i(l ~ "^jr) is the total num- 
ber of flips from state to state 1 inside the observation 
box and, similarly, A/downflips = S^^ (^ ~ '^j,t+i)«it is the 
number of flips from 1 to 0. Of course, trajectories con- 
taining transitions that are not allowed by the dynamical 
rules have Ptraj = and their action is formally infinite; 
we ignore them in what follows. 

We now discuss briefly the extensivity properties of 
the dynamical action. The probability -Ptraj is associ- 
ated with the states of L spins over a time period iobs- 
However, the positions of the spin flip events in time are 
specified with a resolution that depends on the system 
size, TV. Thus, the probability of any trajectory vanishes 
as A^ — > 00 and the action diverges logarithmically with 
N as well as being extensive in L and iobs- We there- 
fore introduce a coarse-graining timescale St. That is, we 
define 



-'traj — / ^ -'tia 

traj'((5t) 



AA„, 



(Ndtr-"'''^' 



.+^fd, 



"'■"-Ptraj, (27) 



where the sum is over trajectories with the same spin 
flips as the original trajectory, but the flips may happen 
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FIG. 8: Distribution of the action in the FA model for L — 60, 
iobs = 320, obtained with N = 180. (Top) Contour plot of the 
joint probability distribution for action density and magneti- 
sation P{ps,m) (obtained from 2 x 10^ independent trajecto- 
ries). The contours are at P{m,pe) = 10^,10,0.1,2 x 10"^. 
The dotted line is the prediction (I34|l . (Bottom) Distribution 
of the action density P{pe) [where pe = f/(Liobs)]. 



within an interval {6t) of their original times. The second, 
approximate, equality indicates that as long as (St) is 
not too large then the probability of all trajectories in 
the sum will be approximately equal, and the number of 
such trajectories is (iV(5t) "'p=. The result is that 



cm 

'traj 



InR 



traj 



= 7Vupfiipsln(l/e ^(50+A/downflipsln(l/(5t) 



-iv-1 



dTY.f^r[il 



nir)e ' + n. 



i=l 



-(boundary terms), 



(28) 



where we have converted the sum to an integral, which 
is valid in the limit of large N . This definition of the 
action is parametrised by the coarse-grained time scale 
{5t). For fixed 5t, distributions oiSl^J- have a peak whose 
position is independent of N and extensive in L and iobs ■ 
Our results do not depend qualitatively on 5t] the nu- 
merical results below are aX 5t = 1. Finally, we have no 
prescription for calculating the boundary terms in the ac- 
tion; these terms are not extensive in Liobs, so we neglect 
them in what follows. 

Results for the distribution of the dynamical action 
density, p£ = £''^*'/(Ltobs), and its joint distribution 



with the box magnetisation, are shown in figure [H We 
see that the distribution of the action density is tightly 
correlated with that of the magnetisation. 

To understand this, we define the density of facilitated 
spins in a trajectory to be 



Nt„ 



nf = {LNt,ksr^ E E/- 



(29) 



Now, facilitated down spins flip with rate e^^ and fa- 
cilitated up spins flip with rate unity. Since tobs 3> i^/c), 
there are many such independent events within the ob- 
servation window, so that fluctuations in the number of 
these events are small. Therefore, 



ATtobs L 



(LTViobs)"' ^ J2frril-' 



T i—1 



T i—1 

-'^upflips — -'^downflips — ^-^^obs^f i 



(l-c)nf,(30) 

(31) 

(32) 



where the approximate equalities indicate that the joint 
distribution of each pair of observables are sharply 
peaked around these values (we have verified this by mea- 
suring joint distributions of the form shown in Fig. [H]). 

Further, each up-spin is associated with two facilitated 
spins. Taking into account spins that are facilitated by 
two separate up spins, we have 



rii 



mti.aj(2- c). 



(33) 



where the approximate equality holds in a similar sense, 
and TOtraj is the magnetisation or activity of the trajec- 
tory of interest, defined in ([5]). 

Hence, taking these results together, we predict that 
the joint distribution of magnetisation and action density 
is sharply peaked around 



P£ = mc{2 - c) [2 + /3 - 2 \n{dt)] , 



(34) 



in accordance with the results of figure [8] We conclude 
that the exponential tail in the magnetisation distribu- 
tion is intrinsically linked with the exponential tail in the 
distribution of the dynamical action. Numerical evidence 
for this result was given in [3[, but without concrete ex- 
planation. 



IV. GENERALISED MODEL 

The definition P^ motivates us to define a new en- 
semble for the dynamics, in which the probabilities of 
trajectories are 



2-lg-f>£.,a 



(35) 



This ensemble has the action distribution 

- lnPb{£) = -[lnPb=i(£) + (6 - !)£]. (36) 
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FIG. 9: (Left) Action distribution in the ensemble with finite b = (1.0, 1.005, 1.001, 1.0015, 1.002, 1.003). The distribution at 
& = 1 is that of Fig. [S] and is shown with symbols. To get data at b > 1 we simply use (|36|l and rescale by a constant for 
convenience (this data is shown as simple lines). (Right) Action distribution with varying R a.t L = 60, tobs = 320, 13 = 1. For 
i? > we use N = 6L = 360 to ensure that data is independent of A^. The behaviour at small R is qualitatively similar to 
the behaviour at small b in that the gradient of the exponential tail decreases; at larger R a secondary minimum appears. The 
inset shows an expanded view of the secondary minimum that is present at R = 0.2. Samples with the action exactly equal to 
zero are omitted from the plot: the probability of this happening is of the order of 1% at _R = 0.2. 



Hence, if Pb=i{£) has an exponential tail with gradient 
a then Pb{£) has a similar exponential tail with gradient 
a — [h ~ 1). Increasing b reduces the gradient, so large 
inactive regions become more common. 

Observables within the 6-ensemble can be measured 
by averaging over the dynamics at 6 = 1 and reweighting 
the resulting trajectories according to their action. See 
the left panel of figure [51 Since the magnetisation and 
action are tightly correlated, increasing h leads to a re- 
duction in the gradient of the exponential tail of P(rn). 
We associate the vanishing of this gradient with the pro- 
liferation of trajectories with large bubbles. In [^, it was 
argued that this proliferation would appear as a first or- 
der transition to a state with large inactive bubbles, if 
an appropriate dynamics could be found to generate the 
6-ensemble. 

For a given L and iobs, a- stochastic process exists 
whose statistics are those of the 6-ensemble (the prob- 
ability of each trajectory is known, and an appropriate 
dynamics can be constructed). This construction is ex- 
tremely laborious, since the transition probabilities de- 
pend explicitly on the size of the space-time box, and on 
the position (and time) within it. Alternatively, one can 
simulate master equations such as those of 9i|, which do 
not conserve probability. Such a simulation requires a 
dynamics in which 'clones' of the system are created and 
destroyed, according to certain rules [25j . 

The connection to experimentally accessible dynamics 
for either alternative is not yet apparent, which limits 
analysis of the predicted first order transition. On the 
other hand, the effects of phase coexistence can be probed 
experimentally by moving towards the second order crit- 
ical point at the end of the phase coexistence boundary. 
As discussed in section IIIIBI above, this critical point is 
approached as T ^ in the FA model. In that limit, the 
order parameter {m) and the surface tension cti vanish 



with the same power of c. The function P{m) scales ac- 
cording to ^ , and it is hard to discern the two coexisting 
phases. 

However, we can modify the FA model by introduc- 
ing an extra process whereby sites can flip from state 1 
to state 0, even when unfacilitated. This new process 
breaks detailed balance, so the resulting model is not ap- 
propriate for an equilibrium physical system, but such 
models are widely studied in a range of physical contexts 
|l2{. The rate for spins to flip from down to up is now 
finite at the critical point, and the new transition is in 
the directed percolation (DP) universality class, where 
the effects of phase coexistence are more apparent. 

To be precise, the generalised FA model has dynamical 
rules 



n^ 


= - 


-^ nj 


= 1 


probability /^e ' , 


rii 


= 1 - 


-^ Tli 


= 


probability /^ -f r(l 



The generalised model has r > since all probabilities 
must be positive; r = is the FA model. For r > then 
the system no longer obeys detailed balance, and /3 can 
no longer be interpreted as an inverse temperature. 

The qualitative behaviour of the model with finite r is 
shown in Fig. 1101 The "dimensionless" parameter that 
determines the effect of r is the product of the death rate 
and the relaxation time of the FA model: 



= 3/3 



i?, 



(37) 



which defines R. At any finite [3 there is a directed per- 
colation transition to an active state, which occurs when 
R is of order unity. This transition is accompanied by a 
diverging static correlation length: 



\ze-d \ri-rj\/£,-DF 



(38) 



where ^dp is the length scale that diverges at the tran- 
sition and z and 9 are DP exponents T2] (the dimen- 
sionality is d = 1). On the other hand, if i? = then 
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FIG. 11: Sample trajectories at _R = 0.1 with conditions 
otherwise similar to Fig. |3] (/3 = 1, L = 60 iobs ~ 320). 
(Left) Sample from centre of distribution. (Right) Sample 
with mtraj ~ {m}/2. Clearly there are more large inactive re- 
gions in these trajectories than in those of Fig. [T] increasing 
R from zero leads to proliferation of large "bubbles" . 



FIG. 10: Sketch of the steady state density in the generalised 
model, as a function of r, for different values of /3 with /3i > 
02 > Pa- The axis r = is the FA model and the axis m — 
is a line of critical points. The dotted line separates the region 
in which the scaling of directed percolation (DP) will apply 
from those in which r can be treated perturbatively [so the 
scaling will be that of the coagulation-diffusion (CD) fixed 
point]. The FA model (r = 0) is the unique case for which 
the critical scaling is coagulation-diffusion; for finite r the 
relevant critical point is DP. 



there is a transition at c = that is controlled by the 
coagulation-diffusion fixed point. In that case we have 



(n„n-jr) - {nirY 



(39) 



for all c. We show in the figure how the crossover into 
the DP critical region is always relevant if R is finite, but 
does not affect the behaviour if i? = 0. We note that 
between two and four dimensions a similar picture holds, 
but the scaling at i? = will be Gaussian 11 1 . 
Turning to the dynamical action, we define 



f< 



(r) ^ 
traj 



M,pflip.ln(e'^)+<tXln(lA) 

dr ^ {/ir[(l - nir)e~^ + n,^(l - r)] + Uirr] 

i 

(40) 

by analogy with ((28| . where we have set the coarse- 
graining time scale 5t — \ for ease of writing. Equa- 
tions (|30ll33p remain true since facilitated spins are still 
equilibrated at a density close to c = {l + e^)~^ (but note 
that the mean density in the system will be less than c: 
the death process reduces the mean density). 

In Fig. [5] we plot the distribution of the action den- 
sity at finite r and compare it with the behaviour of the 
^-ensemble. Since the system does not obey detailed bal- 
ance the TPS procedure becomes inefficient: data at fi- 
nite r is obtained by simple binning and histogramming. 
In the 6-ensemble, the mean action is reduced as the (fi- 
nite) system crosses over from an active to an inactive 
state (this crossover would be a first order phase transi- 
tion in the thermodynamic limit). On the other hand, 
increasing r reduces the gradient of the exponential tail 



and the mean action, as the second order phase transition 
is approached (at Tc). 

The key difference between the FA model and its gen- 
eralisation to finite r is that the surface tension and 
mean magnetisation decouple at finite r. On reducing 
the temperature in the FA model (r = 0, c ^ 0) we have 
tTi ~ (m) ~ c, so the typical bubble size a^ scales in the 
same way as the inverse density of up spins. In contrast, 
as we approach the DP fixed point (/3 finite, r — > Tcrit) 
then we have 



0-1 



{r 



■ Y 



(m) ^ {r - Tcrit) 



/3dp 



(41) 



where (i/dp,/3dp) — (1.1,0.3) are the exponents of the 
DP fixed point [l^] (/3dp should not be confused with the 
parameter (3 that enters the transition probabilities) . As 
criticality is approached, the mean bubble size increases 
much faster than the inverse density since j/dp > /3dp- 
The up spins cluster together, leaving large inactive re- 
gions consisting only of down spins. Introducing the 
death process drives the system towards a phase separa- 
tion into active and inactive regions of space-time. This 
can be seen in the trajectories of Fig. [TT] 

Finally, we note that the critical DP scaling leads to 
a crossover at small R, when the typical size of inactive 
regions exceeds the observation box (criL ~ 1). In that 
case, Pips) acquires a second peak at zero magnetisation 
(see the largest value of R in Fig. ^ . This is analogous 
to the situation for small boxes {cL < 1) in the FA case 
(recall Fig. [5]). In the thermodynamic language, it corre- 
sponds to a second minimum in the free energy, as would 
be observed near a first order phase transition. Here, 
the phase transition is second order, but the finite size 
of the observation box means that the divergence of the 
correlation length is cut off at the size of the observation 
box. 

Thus, we have shown that adding a death process to 
the FA model strengthens the analogy between trajectory 
statistics and the statistical mechanics of phase coexis- 
tence. In the FA model, the only length scale is the in- 
verse density of up spins; if r is finite then a new length 
scale appears: the DP correlation length. This second 
length scale allows the density of active regions to de- 
couple from their separation, leading to richer behaviour 
than that of the FA model (which corresponds to the 
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particular case of the phase coexistence in which the two 
length scales are equal). 

Throughout this article, we have emphasised the broad 
applicability of the idea of space-time phase coexistence 
by focussing on general features of the models, such as 
microscopic reversibility and reducibility of the dynam- 
ics. We expect the behaviour described here to be generic 
in kinetically constrained models with reducible dynam- 
ics 0, QUE 113, Hi]. More generally, the extent to which 
this behaviour can be observed in real structural glass- 
formers remains an open question. 
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APPENDIX A: TRANSITION PATH SAMPLING 

We calculated probability distributions for various ob- 
servables in the FA model: these were obtained by a 
combination of transition path sampling (TPS) [21] and 
umbrella sampling [22] . Here we give a brief description 
of the procedure used. 

Transition path sampling is a Monte Carlo procedure 
applied to trajectories of a dynamical system (a trajec- 
tory is a particular realisation of the dynamics). We 
use it to efficiently sample restricted ensembles of tra- 
jectories. The procedure is well-established [21|, but we 
present some information regarding its application to the 
FA model for completeness. 

Suppose that we wish to sample trajectories with mag- 
netisation m smaller than some reference value mi. The 
simplest way to do this is to generate statistically in- 
dependent trajectories, accepting only those that satisfy 
the restriction m < mi. However, if mi is much smaller 
than the mean of m then this is inefficient. 

Instead, we can generate unbiased trajectories satisfy- 
ing the restriction by deforming a set of initial trajecto- 
ries, as long as the deformations satisfy detailed balance 
with respect to the probability distribution within the 
constrained ensemble of trajectories. Such deformations 
are called TPS moves. One possibility is to take a trajec- 
tory of length iobs and keep only the part with t < (xiobs) 
with < X < 1; we the generate the rest of the trajectory 
using the dynamical rules prescribed for the swtem (we 
use a continuous time Monte Carlo algorithm [29|). If the 
new trajectory has m < m,i then it replaces the old one; 
otherwise the move is rejected and we retain the old tra- 
jectory. This is a "shooting" move [21|. We couple this 



the move with the reverse procedure in which we discard 
the part of the trajectory with t > (xiobs) and regener- 
ate the rest of the trajectory by propagating backwards 
in time (since we have detailed balance then the steady 
state of the FA model is invariant under time-reversal, so 
this is a valid way to generate unbiased trajectories). 

In addition to these moves, we also use "shifting" 
moves [21i] in which we shift the trajectory in time, dis- 
carding the parts of the new trajectory with t < or 
t > iobs- We then regenerate the remaining parts of the 
new trajectory. We find that this combination of moves 
is quite efficient for exploring the restricted ensembles of 
interest. 



APPENDIX B: UMBRELLA SAMPLING 

We use umbrella sampling [22] to calculate probability 
distributions by measuring ratios such as 



P, 



i,i+l — 



P{m < mj+i) 
P{m < mi) 



(Bl) 



where mi and rrii^i arc two cutoffs for the variable m 

with TOi+l < JTli. 

Consider an ensemble of all allowed trajectories for the 
system, with their statistical weights. In order to mea- 
sure -Pi,i+i we sample a restricted ensemble which con- 
tains only trajectories with m < mf, these trajectories 
have the same weights as they would have in the orig- 
inal ensemble. We then measure the probability that 
m < mi+i within the restricted ensemble: this probabil- 
ity is equal to Pij^i. 

Our procedure is as follows: 

1. Start with a representative set of trajectories from 
the unrestricted ensemble. 

2. Choose an ordered set of cutoff magnetisations 
{mi > 7712 > ■ • • > m,p) for which we will calcu- 
late the probabilities Pi^i+i- 

3. Explore the unrestricted ensemble, measuring the 
fraction of trajectories with 777 < 7771. (This is done 
by sampling independent trajectories.) 

4. Once we have a good enough estimate for P(777 < 
7771 ), start a restricted ensemble with trajectories 
satisfying 777 < 7771. Typically we store N^ — 100 
such trajectories. 

5. Explore the restricted ensemble using TPS, mea- 
suring the fraction of trajectories with m < 7772. 
Typically this takes N^ = 100 - 10000 TPS moves 
per ensemble member. 

6. Once we have a good enough estimate for P12, we 
discard all trajectories with 777 > 7772 and replace 
them by trajectories with m < 7772. These trajec- 
tories are generated by continuing the TPS proce- 
dure and accepting all trajectories that satisfy the 
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new constraint. The resulting set of trajectories 
are not statistically independent so we equilibrate 
the new ensemble by allowing it to evolve from the 
biased initial condition. Typically we use around 
Njn TPS moves per ensemble member. We test for 
equilibration by tracking the fraction of trajectories 
with m < 7713, since this will be the quantity that 
we will measure on the next step. 



7. We then repeat steps 5 and 6 for increasingly 
restricted ensembles. At each step, we measure 



P, 



i,i+l- 



Once we have the P(m < mi) and the set of P^^i+i then 
it is simple to reconstruct the probability distribution of 
the observable m. 
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